template< typename T, unsigned int R, unsigned int C >
Matrix<T,R,C>& operator+=( Matrix<T,R,C>& a, const Matrix<T,R,C>& b ) {
  for( unsigned int i = 0; i < R; ++i ) {
    a[i] += b[i];
  }
  return a;
}

template< typename T, unsigned int R, unsigned int C >
Matrix<T,R,C> operator+( const Matrix<T,R,C>& a, const Matrix<T,R,C>& b ) {
  Matrix<T,R,C> temp(a);
  temp += b;
  return temp;
}

template< typename T, unsigned int R, unsigned int C >
Matrix<T,R,C>& operator-=( Matrix<T,R,C>& a, const Matrix<T,R,C>& b ) {
  for( unsigned int i = 0; i < R; ++i ) {
    a[i] -= b[i];
  }
  return a;
}

template< typename T, unsigned int R, unsigned int C >
Matrix<T,R,C> operator-( const Matrix<T,R,C>& a, const Matrix<T,R,C>& b ) {
  Matrix<T,R,C>
    temp(a);
  temp -= b;
  return temp;
}

template< typename T, unsigned int N >
Matrix<T,N,N>& operator*=( Matrix<T,N,N>& a, const Matrix<T,N,N>& b ) {
  return (a = a*b);
}

template< typename T, unsigned int N, unsigned int M, unsigned O >
Matrix<T,N,O> operator*( const Matrix<T,N,M>& a, const Matrix<T,M,O>& b ) {
  Matrix<T,N,O> temp;
  for( unsigned int i = 0; i < O; ++i ) {
    for( unsigned int k = 0; k < N; ++k ) {
      for( unsigned int l = 0; l < M; ++l ) {
        temp[k][i] += a[k][l] * b[l][i];
      }
    }
  }
  return temp;
}

